/**
 * @file BoundaryFunction.cpp
 * @author Lishijie (lsj1018845759@outlook.com)
 * @brief 
 * @version 0.1
 * @date 2021-01-03
 * 
 * @copyright Copyright (c) 2021
 * 
 */
#ifndef _BOUNDARYFUNCTION_
#define _BOUNDARYFUNCTION_
#include "Domain.h"
#include "Point.h"
#include "Mesh.h"
#include <iostream>
#include <vector>
#include <algorithm>
#include <list>
#include <map>
#include <cmath>
#include <functional>

#define TEMPLATE template<Dimension DIM>
enum Boundaryinfo{dirichlet = 1,neumann = 2,robin = 3};
TEMPLATE 
class BoundaryFunction
{
public:
    typedef std::function<double(double*)> Function;
public:
    BoundaryFunction();
    ~BoundaryFunction();
    Boundaryinfo boundaryType(){return _BoundaryType;};
    Function& g(){return _g;};
    std::vector<Boundary<DIM>>& Boundaryvector(){return _Boundaryvector;};
    Function& alpha(){return _alpha;};
    Function& beta(){return _beta;};
protected:
    Boundaryinfo _BoundaryType;
    Function _g;
    Function _beta;
    Function _alpha;
    std::vector<Boundary<DIM>> _Boundaryvector;
};

TEMPLATE 
BoundaryFunction<DIM>::BoundaryFunction()
{
}
TEMPLATE 
BoundaryFunction<DIM>::~BoundaryFunction()
{
}
TEMPLATE 
class Dirichlet:public BoundaryFunction<DIM>
{
    public:
    typedef BoundaryFunction<DIM> _base;
    typedef typename _base::Function Function;
    using _base::_Boundaryvector;
    using _base::_g;
    using _base::_BoundaryType;
    using _base::g;
    using _base::boundaryType;
    public:
    Dirichlet();
    Dirichlet(Boundary<DIM> boundary)
    {
        _g = [](double *p)->double{return 0.0;};
        _BoundaryType = Boundaryinfo::dirichlet;
        _Boundaryvector.push_back(boundary);

    }
    Dirichlet(std::vector<Boundary<DIM>> boundaryvector)
    {
        _g = [](double *p)->double{return 0.0;};
        _BoundaryType = Boundaryinfo::dirichlet;
        _Boundaryvector = boundaryvector;

    }

    Dirichlet(Function g ,Boundary<DIM> boundary)
    {
        _g = g;
        _BoundaryType = Boundaryinfo::dirichlet;
        _Boundaryvector.push_back(boundary);
    }
    Dirichlet(Function g,std::vector<Boundary<DIM>> boundaryvector)
    {
        _g = g;
        _BoundaryType = Boundaryinfo::dirichlet;
        _Boundaryvector = boundaryvector;
    }
    ~Dirichlet(){};
    private:
    Function& alpha(){return _g;};
    Function& beta(){return _g;};
};
TEMPLATE 
class Neumann:public BoundaryFunction<DIM>
{
    public:
    typedef BoundaryFunction<DIM> _base;
    typedef typename _base::Function Function;
    using _base::_Boundaryvector;
    using _base::_g;
    using _base::_beta;
    using _base::_BoundaryType;
    using _base::g;
    using _base::boundaryType;
    public:
    Neumann(){};
    Neumann(Function g ,Boundary<DIM> boundary)
    {
        _g = g;
        _beta = [](double *p)->double{return 1.0;};
        _BoundaryType = Boundaryinfo::neumann;
        _Boundaryvector.push_back(boundary);
    }
    Neumann(Function g ,std::vector<Boundary<DIM>> boundaryvector)
    {
        _g = g;
        _beta = [](double *p)->double{return 1.0;};
        _BoundaryType = Boundaryinfo::neumann;
        _Boundaryvector = boundaryvector; 
    }
    Neumann(Function g , Function beta ,Boundary<DIM> boundary)
    {
        _g = g;
        _beta = beta;
        _BoundaryType = Boundaryinfo::neumann;
        _Boundaryvector.push_back(boundary);
    }
    Neumann(Function g, Function beta , std::vector<Boundary<DIM>> boundaryvector)
    {
        _g = g;
        _beta = beta;
        _BoundaryType = Boundaryinfo::neumann;
        _Boundaryvector = boundaryvector;
    }
    ~Neumann(){};
};

TEMPLATE 
class Robin:public BoundaryFunction<DIM>
{
    public:
    typedef BoundaryFunction<DIM> _base;
    typedef typename _base::Function Function;
    using _base::_Boundaryvector;
    using _base::_g;
    using _base::_beta;
    using _base::_alpha;
    using _base::_BoundaryType;
    using _base::g;
    using _base::boundaryType;
    public:
    Robin();
    Robin(Function g,Boundary<DIM> boundary)
    {
        _g = g;
        _beta = [](double *p)->double{return 1.0;};
        _alpha =  [](double *p)->double{return 1.0;};
        _BoundaryType = Boundaryinfo::robin;
        _Boundaryvector.push_back(boundary);
    };
    Robin(Function g,std::vector<Boundary<DIM>> boundaryvector)
    {
        _g = g;
        _beta = [](double *p)->double{return 1.0;};
        _alpha =  [](double *p)->double{return 1.0;};
        _BoundaryType = Boundaryinfo::robin;
        _Boundaryvector = boundaryvector;
    };
    Robin(Function g , Function alpha ,Function beta ,Boundary<DIM> boundary)
    {
        _g = g;
        _beta = beta;
        _alpha =  alpha;
        _BoundaryType = Boundaryinfo::robin;
        _Boundaryvector.push_back(boundary);
    }
    Robin(Function g, Function alpha ,Function beta , std::vector<Boundary<DIM>> boundaryvector)
    {
        _g = g;
        _beta = beta;
        _alpha =  alpha;
        _BoundaryType =Boundaryinfo::robin;
        _Boundaryvector = boundaryvector;
    }
    ~Robin(){};
};

TEMPLATE 
class BoundaryCondition
{
    public:
    BoundaryCondition(){};
    BoundaryCondition(BoundaryFunction<DIM>* bf){_bflist.push_back(bf);};
    BoundaryCondition(std::vector<BoundaryFunction<DIM>* > bflist){_bflist = bflist;};
    BoundaryCondition(const BoundaryCondition<DIM>& bc)
    {

        for(int i = 0;i < bc._bflist.size();i++)
        {
            BoundaryFunction<DIM>* t;
            Boundaryinfo bfinfo = bc._bflist[i]->boundaryType();
            if(bfinfo == 1)
                t = new Dirichlet<DIM>(bc._bflist[i]->g(),bc._bflist[i]->Boundaryvector());
            else if(bfinfo == 2)
                t = new Neumann<DIM>(bc._bflist[i]->g(),bc._bflist[i]->beta(),bc._bflist[i]->Boundaryvector());
            else if(bfinfo == 3)
                t = new Robin<DIM>(bc._bflist[i]->g(),bc._bflist[i]->alpha(),bc._bflist[i]->beta(),bc._bflist[i]->Boundaryvector());
            else
                std::exit(-1);
            this->_bflist.push_back(t);
        }
    }
    ~BoundaryCondition()
    {
        for(int i = 0;i < _bflist.size();i++)
            delete _bflist[i];
    };
    void add(std::vector<BoundaryFunction<DIM>* > bflist){_bflist = bflist;};
    void add(BoundaryFunction<DIM>* bf){_bflist.push_back(bf);};
    void apply(Mesh<DIM> * m);
    BoundaryCondition<DIM> &operator=(const BoundaryCondition<DIM> & bc)
    {
        for(int i = 0;i < bc._bflist.size();i++)
        {
            BoundaryFunction<DIM>* t;
            Boundaryinfo bfinfo = bc._bflist[i]->boundaryType();
            if(bfinfo == 1)
                t = new Dirichlet<DIM>(bc._bflist[i]->g(),bc._bflist[i]->Boundaryvector());
            else if(bfinfo == 2)
                t = new Neumann<DIM>(bc._bflist[i]->g(),bc._bflist[i]->beta(),bc._bflist[i]->Boundaryvector());
            else if(bfinfo == 3)
                t = new Robin<DIM>(bc._bflist[i]->g(),bc._bflist[i]->alpha(),bc._bflist[i]->beta(),bc._bflist[i]->Boundaryvector());
            else
                std::exit(-1);
            this->_bflist.push_back(t);
        }
        return *this;
    }
    std::vector<BoundaryFunction<DIM>*>& BFlist(){return _bflist;};
    bool is_all_dirichlet()
    {
        bool isalldirichlet = true;
        for(int i = 0;i  < _bflist.size();i++)
        {
            if(_bflist[i]->boundaryType() != Boundaryinfo::dirichlet)
            {
                isalldirichlet = false;
                break;
            }
        }
        return isalldirichlet;
    }
    bool have_dirichlet()
    {
        bool havedirichlet = false;
        for(int i = 0;i  < _bflist.size();i++)
        {
            if(_bflist[i]->boundaryType() == Boundaryinfo::dirichlet)
            {
                havedirichlet = true;
                break;
            }
        }
        return havedirichlet;
    }
    private:
    std::vector<BoundaryFunction<DIM>*> _bflist; 

};
template<>
void BoundaryCondition<2>::apply(Mesh<2>* m)
{
    bool isalldirichlet = true;
    for(int i = 0;i  < _bflist.size();i++)
    {
        if(_bflist[i]->boundaryType() != Boundaryinfo::dirichlet)
        {
            isalldirichlet = false;
            break;
        }
    }
    if(isalldirichlet)
        return;
    else
    {
        for(int k = 0;k < m->_BndEdge.size();k++)
        {
            Dofs<2> dof1 = m->DofsofIndex(m -> _BndEdge[k][2]);
            Dofs<2> dof2 = m->DofsofIndex(m -> _BndEdge[k][3]);
        for(int i = 0;i < _bflist.size();i++)
        {
            Boundaryinfo BndType = _bflist[i]->boundaryType();
            if(BndType != Boundaryinfo::dirichlet)
            {
                std::vector<Boundary<2>> BF = _bflist[i]->Boundaryvector();
                for(int j = 0;j < BF.size();j++)
                {
                    Boundary<2> B = BF[j];
                    if(B.is_on_boundary(dof1) && B.is_on_boundary(dof2))
                    {
                        if(BndType == Boundaryinfo::neumann)
                            m->_BndEdge[k][0] = -2;
                        else if(BndType == Boundaryinfo::robin)
                            m->_BndEdge[k][0] = -3;
                        else
                            std::exit(-1);
                    }
                }
            }
        }
        }
    return;
    }
};
#undef TEMPLATE 
#else
#endif

